data = readtable('../../data/main_VAR/UK_source_pre1946.xlsx', 'Sheet', 'Full_Data_for_VAR');
date = data.Year;

startdatenum = 1729;
enddatenum = 1946;
middatenum = 1794;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
middate = find(date == middatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.Year >= startdatenum - 1), :);
data = data((data.Year <= enddatenum), :);

taxrevgdp = data.revtogdp(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.revtogdp(1); % government tax revenue to GDP ratio
spendgdp = data.spendingtogdp(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.spendingtogdp(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.realGdpGrowth(2:end); % real gdp growth (in logs)
inflation = data.inflation(2:end); % growth in GDP price deflator (in logs)
ynom1 = data.short_r(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.rate_10year(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

dy = data.dy(2:end);
pdm = log(1 ./ dy); % log price/dividend ratio

%% convenience yield

cy_short1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck5:ck45'); % 1873-1913
cy_short2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck57:ck63'); % 1925-1931

cy_short = [mean(cy_short1) * ones(1873 - startdatenum + 1, 1); cy_short1; mean(cy_short2) * ones(1925 - 1914, 1); ...
                    cy_short2; mean(cy_short2) * ones(enddatenum - 1931, 1)] / 100;

cy_long1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck5:ck45'); % 1873-1913
cy_long2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck57:ck63'); % 1925-1931

cy_long = [mean(cy_long1) * ones(1873 - startdatenum + 1, 1); cy_long1; mean(cy_long2) * ones(1925 - 1914, 1); ...
                   cy_long2; mean(cy_long2) * ones(enddatenum - 1931, 1)] / 100;

% cy is zero before 1794;
cy_long = [zeros(middate - startdate + 2, 1); cy_long(3 + middate - startdate:end)];
cy = cy_long;

ynom1 = ynom1 + [zeros(middate - startdate + 1, 1); cy_short(3 + middate - startdate:end)];
ynom10 = ynom10 +cy_long(2:end);

ynom1 = log(ynom1 +1);
ynom10 = log(ynom10 +1);

divgrm = diff(log(data.dy) + data.stockindex);

gdebt = data.debttogdplevel;

% calibrate cy/gdp ratio mean
cy_data = mean([(cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2);
                (cy_long2) / 100 .* gdebt(find(date == 1925) + 1:find(date == 1931) + 1)]);

datew0 = (1728:1946)';
seig_rev = zeros(size(datew0));
seig_rev(find(datew0 == 1795):find(datew0 == 1946)) = cy_data;

taxrevgdp = taxrevgdp + seig_rev(2:end);
taxrevgdp0 = taxrevgdp0 + seig_rev(1); % use gdebt in year 1729 as the value in 1728;

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1; % spread between 10-yr and 1-yr yields
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% log div/gdp ratio growth
divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



tts(1) = taxrevgdp(1);
gts(1) = spendgdp(1);

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

%% 
run ../../tools/run_var_estimate.m

save(['MAT/res_annual_VAR_nodebt_cy_robust_partial.mat'],'-regexp', '^(?!(option|Globaloption)$).');
